Agglomeration of Oppositely Charged Particles in Nonpolar Liquids 
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We study the aggregation of insulating electrically charged spheres suspended in a nonpolar 
liquid. Regarding the van der Waals interaction as an irreversible sticking force, we are especially 
interested in the charge distribution after aggregation. Solving the special case of two oppositely 
charged particles exactly, it is shown that the surface charges either recombine or form a residual 
dipole, depending on the initial condition. The theoretical findings are compared with numerical 
results from Monte Carlo simulations. 

PACS numbers: 83.10.Rs, 05.10.Gg, 81.07.Wx, 45.70.-n 



I. INTRODUCTION 



Fine powders with particles on the micrometer scale 
play an increasing role in diverse technological appli- 
cations, raiiging from solvent-free varnishing to inhal- 
able drugs A major problem in this context is the 
tendency of the particles to clump due to mutual van 
der Waals forces leading to the formation of aggre- 
gates. In many applications, however, the aggregates 
should be sufficiently small with a well-defined size dis- 
tribution. A promising approach to avoid clumping is 
to coat the powder by nanoparticles. The small par- 
ticles act as spacers between the grains, reducing the 
mutual van-der- Waals forces and thereby increasing the 
flowability of the powder. However, the fabrication of 
coated powders is a technically challenging task since the 
nanoparticles themselves have an even stronger tendency 
to clump, forming large aggregates before they are de- 
posited on the surface of the grain. One possibility to 
delay or even prevent aggregation is the controlled use 
of electrostatic forces. As shown in Ref. T] this can be 
done by charging the nanoparticles and the grains oppo- 
sitely. On the one hand the repulsive interaction between 
equally charged nanoparticles suppresses further aggre- 
gation once the Coulomb barrier between the flakes has 
reached the thermal energy On the other hand, at- 
tractive forces between the nanoparticles and the grains 
support the coating process. 

The coating process is most easily carried out if both 
fractions of particles are suspended in a liquid (see 
Fig. This type of coating processes requires the use 
of a nonpolar liquid such as liquid nitrogen. In contrast 
to colloidal suspensions in polar liquids, the charged par- 
ticles suspended in liquid nitrogen are not screened by 
electrostatic double-layers. Both the large and the small 
particles are insulators so that the charges reside on their 
surface. By choosing different materials and charging 
them triboelectrically, it is possible to charge the two 
particle fractions oppositely in a single process f^. 

What are the morphological properties of the coated 
surface? Do the nanoparticles reach their countercharges 





FIG. 1: Schematic experimental setup of the coating process: 
Large and small particles are suspended in a nonpolar liquid. 
Charging them oppositely the small particles are preferentially de- 
posited at the surface of the large particles. 
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FIG. 2: If two oppositely charged particles aggregate due to van 
der Waals forces, the charges may either recombine or form a resid- 
ual dipole. 



exactly or do they attach elsewhere on the surface, form- 
ing residual dipoles, as sketched in Fig. In order to 
address these questions we consider a simplified situa- 
tion, where pointlike particles are deposited on the sur- 
face of a plane, representing the surface of an infinitely 
large spherical particle (see Fig.l^J. One or several posi- 
tive charges are located on the planar surface of the big 
particle, attracting negatively charged pointlike particles 
inserted far away. For simplicity we assume that the par- 
ticles are inserted one after another so that mutual inter- 
actions during the deposition process can be neglected. 
Similarly, we assume that the hydrodynamic interactions 
between the particle and the plane can be ignored. Thus 
the small particle is subjected to Coulomb forces and 
Stokes friction as well as Brownian motion. As shown 
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FIG. 3: Simplified model for electrostatically supported coating: 
A negatively charged pointlike particle is attracted by a positive 
charge located at the planar surface of an infinitely large spherical 
particle. 



in Ref . ^ , the damping time of suspended nanoparticles 
is so short that on the time scale relevant for the coat- 
ing process their motion can be assumed as overdamped, 
i.e., inertia can be neglected. The van der Waals interac- 
tion is interpreted as a purely adhesive force, i.e., once a 
particle touches the surface of the large particle it sticks 
irreversibly. 

In this study we show that a certain fraction of the 
particles exactly reach and compensate their counter- 
charges. The remaining particles are distributed around 
the charges, partly decaying with distance, partly as a 
constant background. The fraction that exactly reaches 
the countercharges is determined by the interplay be- 
tween the magnitude of the charges, the density of the 
charges, and the diffusion constant. It can be used as 
a measure to what extent a predefined structure of pos- 
itive charges at the surface survives during the coating 
process. 



II. THEORETICAL PREDICTIONS 



A. Formulation of the problem 



In order to study the deposition process analytically, 
one has to solve the equation of motion of the parti- 
cle subjected to Coulomb forces, Stokes friction, as well 
as Brownian motion. We work in the overdamped limit 
which can be justified as follows. On the one hand, the 
viscous motion of the particle decays exponentially on a 
time scale 



^rclax 



(1) 



where 77 is the viscosity of the fluid, a is the particle 
radius, and m the particle mass. On the other hand, the 
typical time scale for a particle to diffuse thermally by 
its own diameter is given by 



^diff 



(2) 



As shown in Ref. under typical experimental condi- 
tions in liquid nitrogen treiax is always much smaller than 
^diff, even if the particles are as small as 1 nm. There- 
fore, on time scales larger than tdiff, the particle performs 
a random walk guided by the balance of Coulomb forces 
and Stokes friction. Such a motion can be described by 
the Langevin equation 



d 



Fcjr) 



(3) 



where r is the position of the particle, Fc{r) is the 
Coulomb force acting on it, and ^(i) is a white Gaus- 
sian noise with the correlations 



5^j5{t~t'). 



(4) 



Equivalently, one may formulate the problem in terms of 
a Fokker-Planck (FP) equation which describes the 
temporal evolution of the probability distribution P(r, t) 
to find the particle at point r and time t. The FP equa- 
tion has the form 



-V-j(r,i), (5) 
where 

J (r, t) = -DVPif, t) + v{f)P{r, t) (6) 
is the probability current, D the diffusion constant, and 



v{r) = 



9 



Q r 



is the particle velocity in the overdamped limit. 
Rescaling space and time by 



(7) 



D 



r , 



t 



t 



(8) 



and suppressing the arguments r , t we obtain the 
parameter-free dimensionless equation 



d_ 

dt 



P ^ V^P - VP • u - P (V • w) 



(9) 



where u 



-rjr' 



B. Solution of the Fokker-Planck equation 



In what follows we consider a pointlike particle inserted 
at a finite distance z—L from the plane with random 
coordinates x and y, as shown in Fig. ^ The particle 
diffuses guided by the Coulomb force until it touches the 
wall at 2;=0, where it sticks irreversibly. Our aim is to 
compute the probability distribution p(r) for the particle 
to touch the wall at a distance r=\/ -\- ip-. To this end 
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FIG. 4: A thermally diffusing particle inserted randomly at finite 
distance z=L is attracted by the target charge at the origin. Once 
the particle touches the target plane 2=0, which represents the 
large particle to be coated, it sticks to its surface due to van-der- 
Waals forces. The single-particle problem can be recast as a quasi- 
stationary flow of a large number of particles inserted continously 
along the source plane which are removed whenever they touch the 
target plane. 

we consider the problem as a quasi- stationary process, 
where many particles are continuously introduced at the 
source plane z=L and removed at the target plane z=0 
(see Fig. 0J| . Thus the probability distribution Pg (r) to 
find a particle at position r is a solution of the stationary 
FP equation 

V^F, -VP, •u-P,(V-?2) = (10) 
together with the boundary condition 

^«(^1L=o = o (11) 

and an appropriate source term at z=L. Taking the limit 
L ^ oo this problem can be solved exactly in two (2D) 
and three (3D) spatial dimensions (see Appendices IXIand 
ig. In the original variables the stationary probablity 
distribution is given by 

^'^"^ " { z{l + l/2r) in 3D . ^^^^ 
The probabilitiy current in 2D is given by 

/ = ^^/^' = ( -^r~iyr^ ) ' (13) 

while in 3D one arrives at a more complicated expression 

--(-i)(;)-^;^(»)' <») 

Remarkably, in both cases the vector field j{r) exhibits 
a separatrix between a region, where the fiux lines reach 
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FIG. 5: Flux lines of the vector field j{r) in the x—z plane for 
the three-dimensional problem. The separatrix between flux lines 
terminating at the target charge and those terminating elsewhere 
at the wall is represented as a bold line. In the two-dimensional 
case one obtains a qualitatively similar flow field. 



the target charge at the origin, and another region, where 
they terminate elsewhere on the surface (see Fig.lSJ. As 
shown in Appendix [0 it is even possible to calculate the 
separatrix exactly. Note that the fiux lines and the actual 
trajectories of the particles have a different meaning. 

The density of particles p{r) reaching the wall at dis- 
tance r from the target charge is proportional to the 
normal component of the fiux jz{^\z=o- shown in 
Appendix^ we obtain 

[1 + 25(0 in 2D 

P^"^' " \ 1 + l/2r + f <5(r) in 3D . 

Concerning the problem of several target charges, we note 
that in 2D the problem is still analytically solvable. For 
a given charge density A(r) on the target line, the density 
of charges reaching the wall becomes 

p(r) = l + 2A(r) . (16) 

This result stems from the fact that in 2D the station- 
ary probability distribution Ps{r)—z, Eg. (fT^ . is inde- 
pendent of the charge position, or more generally, inde- 
pendent of the distribution A. This can be verified by ex- 
plicitly inserting the solution in Eq. (|10(l . The Laplacian 
V^Ps vanishes so that the resulting equation becomes lin- 
ear in the Coulomb potential that is contained in u. The 
proposed solution for P, solves the equation for any point 
charge, and due to this linearity, it solves the equation for 
any charge distribution. In the physically more relevant 
three-dimensional case, however, we resort to numerical 
methods to study the effect of more charges. 
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FIG. 6: Boundary conditions in the numerical simulation: The 
picture shows a two-dimensional cut through the three-dimensional 
simulation setup. The x-axis points perpendicular to the shown 
plane towards the reader. 



III. NUMERICAL RESULTS 



A. Implementation 



The discretized equation of motion for a single particle 
can be derived from the Langevin equation and is 
given by 



df 



-^-dt + V2Ddtn 



(17) 



where Q = Dq^ /{STreeokBT) is a constant containing 
the diffusion constant D, dielectric constant eeo, parti- 
cle charge q, and thermal energy fcsT. TZ is a vector 
of Gaussian distributed random numbers with unit vari- 
ance, representing Brownian motion of the particle. As 
in the previous section, this equation can be made di- 
mensionless by rescaling space and time by r — s- -^f* and 

t j^t, leading to 



(18) 



As in the FP equation ®, no free parameter is left in 
this equation. Thus, apart from the units of space and 
time, the solution is universal. The numerical integration 
of Eq. (|18|l can be performed easily for a large number 
of representations of the Brownian motion TZ and initial 
positions at the source plane on a workstation. 

A two-dimensional cut through the three-dimensional 
simulation setup is shown in FigEl One or several charges 
are fixed at the y-axis. Small particles start diffusing 
from the plane (—50 < x,y < +50, z — -1-5) towards the 
y-axis. In order to avoid diffusion of particles too far 
away from the interesting region, the simulation volume 
is confined by walls at x,y = ±70 and z = +70. Particles 
touching this walls stick irreversibly. 



Comparision of numerical results and 
theoretical predictions 



First we want to show that in case of a single charge 
fixed at the wall our numerical simulations reproduce 
the exact analytical solution. This has to be done, since 
the boundary conditions of our simulation setup (diffus- 
ing particles start from a plane with finite extent, not 
far away from the absorbing wall; additional absorbing 
walls confine the simulation space) are different from the 
boundary conditions of the Fokker-Planck-equation seen 
above. 

By running the simulation without any attracting 
charges fixed at the walls (i.e. particles only diffuse) one 
can compare simulation results to the homogeneous dis- 
tribution of particles on the wall one would expect from 
solving the FP equation without any attracting charge. 
The numerical results show a homogeneous distribution 
of particles hitting the wall in a sufficiently big region. 
However, approximately 14% of the particles hit the ad- 
ditional walls surrounding the simulation space, leading 
to a reduced influx of particles to the z=0 plane. This 
will be taken into account in the following graphs by am- 
plification of the incoming particle fluxes to compensate 
for this loss of particles. 

In order to check the numerical results for the case of 
one attracting charge fixed at the wall, we calculate the 
influx of particles in a circular region of radius R around 
the attracting charge. From equation (|15|l . the influx is 
given by 



1 \ TT 

— 27rrdr + -5{r) = -kR^ 
2r 2 ^ ' 



ttR 



Fig.El shows the influx obtained in the simulation for dif- 
ferent bin sizes R compared to the theoretical influx after 
subtraction of the homogeneous background influx. As 
one can see, both are in excellent agreement. 



C. Numerical results for serveral target charges 



What changes, if serveral attracting charges are fixed 
at the wall? Fig.|Slshows simulation results for up to four 
charges fixed at the wall. The dimensionless Langevin 
equation is now given by 



dr^Y.- 



[r - 



(20) 



where n is the number of charges. The charges are always 
located on the y-axis, separated by a distance of two 
dimensionsal units. The boundary conditions are still 
the same as shown in Fig. El Fig. |H1 shows the density 
distribution of incoming particles in a small strip around 
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FIG. 7: Influx of particles onto the wall within a circular shaped 
region around the fixed charge. The straight line represents the 
analytical solution: ^+itR. In both cases, the constant background 
is subtracted. 



the y-axis {\y\ < 0.5). In each graph the dark Hne shows 
simulation results, while the light curve is given by 



p{x = 0,y,z = 0) = I 



1 



^ 2|y 



(21) 



where n is the number of charges fixed at the y-axis. 
As one can see, the assumed superposition of the 1/r- 
shoulders from the single charge solution fits suprisingly 
well, since this kind of superposition is not a solution of 
the FP-equation. However, one can see that in the case 
of three and four charges, the density of particles near 
the outer charges is slightly overestimated. 

The relative strength of (5-peaks for different numbers 
of charges is shown in Fig.|^ Again, this data is obtained 
by placing a circular bin around each charge location and 
extrapolating the strength of the delta peak from a van- 
ishing bin radius. As one can see, the amplitude of the 
i5-peak grows by increasing the number of charges. Also, 
charges in the center of line always collect more particles 
than charges on the edge. 



D. Agglomeration of equally sized spherical 
particles 



So far the limit of one infinitely large particle and one 
small particle was examined in detail. Now we address 
the other limit, namely two spherical particles of the 
same size. 

The rotational degree of freedom must be taken into 
account when the two particles are of equal size. As 
for the translational degree of freedom, we assume that 




FIG. 8: Density profile of particle influx for one (a), two (b), 
three (c) and four (d) charges fixed on the y-axis. The bold curves 
show data from our simulations, the ligt curves are computed as 
described in the text. 
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FIG. 9: Relative strenght of the delta peaks for one (1), two (2), 
three (3/ for the inner, 3// for the outer peak) and four (4/ for the 
inner, 4// for the outer peak) charges. 



the overdamped limit is valid also for the rotational mo- 
tion. Again, the rotational motion can be described by a 
Langevin equation, i.e.. 



d 

at' 



3M 



(22) 



where (j) describes the position of the charge relative to 
the z-axis, M is the excerted torque on the particle due 
to Coulomb forces, and ^ means the brownian rotational 
displacement. 

When the system consists of two particles, each carry- 
ing one charge (of opposite sign), that are located at a 
specified point on the surface, then the Coulomb force 
gives a contribution to the rotation as well as to the 
translation of the particles. The Coulomb force tends 
to rotate the two particles so that the charges approach 
the common axis of the particles, the minimum distance. 
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from a regime where it approaches the maximal value. 
This allows one to estimate the typical dipole length, in 
experiments as well as in simulations. 

In reality, particles are not perfectly circular or spher- 
ical. The roughness of the surface of a particle has the 
qualitative implication that rotations become increas- 
ingly difficult with increasing roughness. As a conse- 
quence, we expect a shift of the cross-over to lower ther- 
mal energies. In this sense, the oberserved point can 
be considered as an upper limit for the range where the 
Coulomb forces dominate. 



FIG. 10: Agglomeration of two spherical particles: Average dipole 
strength in units of particle radii a as function of the particle radius 
in units of the Bjerrum length Ig. For particle radii bigger than 
the Bjerrum length, the position of charges during agglomeration 
is randomized due to Brownian rotation. 

Rotational Brownian motion tends to randomize the 
orientation of the particles. This means that it is the rela- 
tive strength of Coulomb force and Brownian motion that 
decide whether the two charges find each other upon the 
collision, or whether a permanent dipole pertains. The 
natural scale for Coulomb energy is the energy of the two 
charges two particle radii apart Ec = /2a. Rotational 
Brownian motion is controlled by thermal energy, and if 
ksT » Ec the orientation of the particles is completely 
random when they collide. The average distance of the 
resulting dipole can in this case, by a simple numerical 
integration over the surfaces of the disks or the spheres, 
be calculated. The numerical values are: 2.26a (2D) and 
2.33a (3D), i.e., slightly more than two particle radii. In 
the other limit of vanishing Brownian motion the aver- 
age distance is of course zero. Further, a crossover in 
the average distance of the dipole is expected roughly 
at ksT/Ec = 1. In colloidal sciences, the crossover be- 
tween regions dominated by coulomb interaction on the 
one hand and by thermal diffusion of particles on the 
other hand is often determined by the Bjerrum length 

Ib . (23) 

47reeofc_Br 

On distances smaller than the Bjerrum length, interac- 
tion of particles is guided by coulomb interaction, while 
on larger distances diffusion dominates. If all colloidal 
particles carry identical charges, the suspension is sta- 
bilized by Coulomb repulsion if the particle radii are 
smaller than the Bjerrum length. 

The values of the limits and the location of the cross- 
over are verified by numerical experiments as shown in 
Fig. IIUI Experiments are performed for differently sized 
particles and for different charges in 2D. The data sets 
collapse since only the ratio of the Coulomb energy and 
the thermal energy is relevant. This universal curve 
shows a cross-over, seperating a regime where the resid- 
ual dipole moment increases lineary with thermal energy. 



IV. CONCLUSIONS 

The use of electrostatic forces strongly affects the dy- 
namics of agglomeration processes, for instance coating of 
micrometer particles with nanometer particles. In par- 
ticular, it is still an open question to what extent the 
exact locations of charges on the particles have impact 
on the resulting agglomerated structures. In this context 
the subquestion arises whether recombination of two par- 
ticles carrying opposite charge gives a pertaining dipole 
or whether the charges cancel out. To simulate large sys- 
tems of suspended charged particles, knowledge about 
charge recombination allows for proper approximations, 
which in turn may increase the range of many-particle 
simulations. 

A detailed study of two idealized situations have been 
presented. First the case of one very large particle with 
a target charge, which attracts a much smaller particle 
with one oppsite charge. Formulating the problem in 
terms of a Fokker-Planck equation, the stationary solu- 
tion is found analytically in 2D and 3D. The findings are 
supported by numerical simulations, and the key results 
are as follows. Interpreted as a hitting probability distri- 
bution, there is a finite fraction of charged particles that 
exactly recombine with the target charge, giving rise to 
a (5-peak in the distribution. This is valid in 2D as well 
as in 3D, although the strength of the (S-peak is some- 
what larger in 2D. However, in 3D there is an additional 
contribution decaying like l/2r, where r is the distance 
from the target charge. This is the main difference be- 
tween 2D and 3D, and it is the reason why superposition 
of the solution is possible in 2D, but not in 3D. Further, 
there is a constant background term in the probability 
distribution whose physical reason lies in the diffusion of 
the particles. 

Numerics is performed to test the analytic results, and 
to provide qualitative insight into the situation of hav- 
ing more target charges in 3D, which from the point of 
view of applications is very important. In this case the 
superposition of a homogenous background plus a l/2r 
shoulder located around each fixed charge fits the numer- 
ical data suprisingly well. Deviations from the superposi- 
tion of single charge solutions are visible in an increasing 
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strength of delta peaks. Thus, an increased number of 
fixed charges increases the probabihty for a compensation 
of fixed charges and incoming particles. 

The second idealized situation is the study of two par- 
ticles of equal size recombining. In this case the rota- 
tional degree of freedom is included. From physical rea- 
soning and from simulation, we demonstrate that there 
is a cross-over when Coulomb energy equal to thermal 
energy. Larger thermal energies means a regime where 
dipoles are created with random distance. For lower ther- 
mal energies the average length of the dipole goes to zero 
in the limit, but for finite ratios of the energies, the aver- 
age length increases rapidly with the energy ratio. This 
implies that for simulations or experiments in this range, 
the detailed localization of charges on the particles must 
be taken into account for a correct description of the 
physics. 
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APPENDIX A: SOLUTION OF THE 
FOKKER-PLANCK EQUATION IN TWO 
DIMENSIONS 



Inserting these expressions, the stationary Fokker-Planck 
equation with u — f/r^ reads 



92 d 



d 



+ r— + {r^+r)—-l]P{r,ip)^0. (A5) 



Using the Ansatz 

one obtains two separate equations 



(A6) 

+ C)QM = (A7) 



& 



d 



ITT + + - (1 + Cr) R{r) - 0, (A8) 



dr 



where C is the common eigenvalue. Symmetry requires 
that Q{uj) is an even function, and the possible solutions 
of the angular equation (jA7|l are given by 



Q{^) — cos{llhp) with C = uj^ 



(A9) 



As the absorbing wall at z—Q imposes the boundary con- 
dition Q(-7r/2) = Q(+7r/2) = we have 

= 1, 3, 5, . . . . 

Inserting C = uj"^ into equation ljA8|l we find the solutions 

uj = l: R(r) = Ar + B{r - 1)6^/'^ 
oj = 3:R{r) = A(20r^ + 8r^ + r) 

+B{60r^ - 36r2 + 9r- l)e^^'' 



In order to solve the stationary Fokker-Planck equation 



(Al) 



with u= — r/r^, we consider a stationary flow of par- 
ticles inserted at infinity and removed whenever they 
touch the wall (cf. Sec. IIIBjl . Let us first consider 
the two-dimensional case in the x, 2-plane, where the z- 
axis denotes the direction perpendicular to the surface 
of the grain (see Fig. 0)) . Introducing polar coordinates 
z=r cos if and x=r sin tp the gradient acting on a scalar 
/(r, (p) and the divergence of a vector field u{r, (p) are 
given by 



V/ = er 



dl 
dr 



idi 

r dip 



^ ^ d Id 1 

V-U = T^Ur + -^U^ + -Ur , 

or r oip r 



(A2) 
(A3) 



where the last term is due to the covariant derivative 
of vector fields in polar coordinates. In particular, the 
Laplacian is given by 



52 1 92 Id 



Qj,2 j,2 Q^p2 J, Qj, 



(A4) 



In general these solutions are of the form 

Ahir)+Bhir)e'/\ (AlO) 

where /i and /2 are polynomials of degree lu. 

Far away from the target charge the probability dis- 
tribution P(r, (fi) will not be influenced by the Coulomb 
force. Since it is assumed that particles are inserted ho- 
mogenously at large distance, we therefore expect a linear 
asymptotic behavior 



P(x, z) cx z for z ~* 00. 



(All) 



Obviously this condition can only be satisfied for lu = I. 
Furthermore, since P must be positive for small r, we 
find B=0. Hence the solution of the two-dimensional 
Fokker-Planck equation simply reads 



P{x,z) ^ A - z 



(A12) 



independent of the surface charge, where A is a normal- 
ization factor equal to the rate of inserted particles per 
unit area at the source plane. Setting A=l the corre- 
sponding probability current j = Pu — VP for r>0 is 
given by 

(A13) 



3 = ex - 1 
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APPENDIX B: SOLUTION OF THE 
FOKKER-PLANCK EQUATION IN THREE 
DIMENSIONS 



APPENDIX C: FLUXLINES OF THE 
PROBABILITY CURRENT j 



Following the previous calculation, we solve the 
Fokker-Planck equation 1A1|I by first transforming it to 
spherical coordinates 



X = r sin cos ip 
y = r sin sin ip 



z 



r cos . 



In these coordinates the gradient acting on a scalar is 
given by 



df 



Idf 



1 df 



V/ = e,. • — + e^i • -— + 

or r ov r sm v 0(p 



while the Laplacian takes the form 



92 2d 1 92 



1 



d 



1 



(Bl) 



92 



Using again a separation ansatz 

we are led to the equations 
/ 92 1 9 1 92 



' tan d dd sin d dtp'^ 
(B2) 

(B3) 



V9i92 taiidd^ sin2^9(^2 



+ C Q(^,^) = 0(B4) 



,9^ 

9r2 



{2r+l)^-C]R{r) = (B5) 



As illustrated in Fig. 0] the system is invariant under ro- 
tations around the z-axis. Thus, the solution will only 
depend on r and hence Q{'d, ^) = Q{'&)- 

Solving the angular equation, the general solution can 
be expressed in terms of Legendre polynomials. However, 
for large r we expect the solution to be independent of the 
Coulomb field, i.e., linear in z. Therefore, the only solu- 
tion of the angular equation, which satisfies the boundary 
condition Q{tt/2) = 0, turns out to be Q^d) — cos-Q with 
the eigenvalue C = 2. The corresponding radial equation 
has the solution 



R{r) = A{2r + 1) + B{2r - 1) e 



l/r 



(B6) 



Since P(r, ??, ip) has to be non-negative for small r the 
second term has to vanish, i.e., B — Q. Choosing A — 1/2 
the physically meaningful solution reads: 

P=(l + l)r.cos^=(^l + i-)z. (B7) 

The corresponding probability current j = Pu — \7 P for 
r > is given by 



1 \ ^ (1 



(B8) 



In two dimensions the trajectories of the vector 
field l|A13|l can be obtained by solving the differential 
equation 



dz 
dx 

leading to the solution 



z{x) = ±' 



(CI) 



(C2) 



where c is an integration constant labeling different 
curves. For the separatrix the slope at the origin 



dz(x) 



dx 



= ± 



(C3) 



vanishes, i.e. c = ±1. Selecting the physically meaningful 
branch the separatrix is given by 



z{x) = ""^f^^ (0 < X < 1) . (C4) 

In three dimensions the separatrix can be calculated in 
the same way. Because of rotational invariance in the xy- 
plane we set y = so that the trajectories of the vector 
field (|B8|) obey the differential equation 



where = 



dz 
dr 

The solution reads 
z{r) 



dz _ 2r^ + + (1 + r)z'^ 
dx (1 -t- r)xz 

' + z2, or equivalently 

2^4 ^j,3 ^ +r)z2 



r.z(2r2 + 2r+l) 



C + 2r + 2r2 
1 + 2r + 2r2 



(C5) 



(C6) 



(C7) 



where C is an integration constant. Since z'(0)=0 implies 
C—0 the separatrix in original coordinates is given by 



z{x) = ^2{x^-x^) + V^ 



(C8) 



APPENDIX D: PARTICLE DENSITY AT THE 
WALL 



The main quantity of interest, which can be calculated 
from the probability density flux, is the distribution of 
the particles that hit the wall. Due to the rotational 
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symmetry of the configuration around the charged parti- 
cle, this density, p{r), is a function of the radial distance r 
from the charge only. The density is equal to the normal 
component of the flux at the wall 



(Dl) 



For all points on the wall, except for r=0, it follows by 
direct insertion into Eqs. I|A13(I and (|B8|I . that 



p{r) 



in 2D 

l/2r in 3D 



(r > 0) 



(D2) 



A certain fraction of the particles hit the target charge 
directly at r—Q^ thus giving a (5-peak contribution at this 
point. In 2D the strength of the peak is 



tt/2 



lim 



Similarily in 3D, by taking the unit normal vector as 
fi — (sin •& cos Lp, sin § sin 93, cos ip) , one obtains 



fi ■ j R dip 



-7r/2 



n/2 



lim I — 



lim — 27r 
2 ' 



d-d / dipR^ sin-dn ■ j 



J 



+ R 



7r/2 



COS d sin § d-d 



(D4) 




where the integration is taken over a half-sphere around 
the origin with radius R. Combining these results we 
arrive at 



p{r) 



lim 21 R 



1 + 2S{r) in 2D 

1 + l/2r+ f(5(r) in 3D 



(D5) 



where n—Cz • cos p + Cx ■ sin ip is the normal vector on 
a half-sphere over which the integral in-flux of particles 
is calculated. The actual in-flux at r—0 is found in the 
limit r^O. 



In the 2D solution it is understood that the ^-function 
integrated over the target line gives unity. Likewise in 
3D, integration of the (5- function over the target plane 
gives unity. 
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